
/* This file contains functions which implement a sequential algorithm to model the steady state equations for a metbolic network */

// includes, system
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "steadyState.h"

float max(float* input, int size);
//void ReadFile(char* filename, void* input);

int main(int argc, char** argv) {

//Declare Variables
//Inputs
float epsilon = 0.01f;
float fwdRate[NUM_EDGES] = {1,1,1,1,1,1,1,1};      // forward rate constant - kf
float revRate[NUM_EDGES] = {0,1,1,1,1,0,1,0};      // reverse rate constant - kr
float metaCnc[NUM_NODES] = {1,1,1,1,1,1};      // metabolite concentration - m
float fullStoic[NUM_EDGES * NUM_NODES] = { 1,-1,-1,0,0,0,0,0,
				           0,0,2,-1,-1,0,0,0,
				           0,1,0,0,-1,0,1,0,
				           0,1,0,0,0,0,-1,-1,
				           0,0,0,0,1,-1,0,0,
				           0,0,0,1,0,-1,0,0};  // full stoichiometric matrix - S
float tolerance = 0.00001f;
unsigned int timeStep;
float t = .001f;
//int numReactions = ;
//int numMetabolites = ;
int numReactions = NUM_EDGES;
int numMetabolites = NUM_NODES;

//Outputs
float scaleFactor[numMetabolites];  // volume scaling factor - a
float fwdFlux[numReactions];      // forward flux - vf
float revFlux[numReactions];      // reverse flux - vr
float totFlux[numReactions];      // total flux - vt
float metaCncDiff[numMetabolites];  // change in metabolite concentration - dm

float fwdStoic[NUM_EDGES * NUM_NODES];  //substrate stoichiometric values in forward direction - Sf
float revStoic[NUM_EDGES * NUM_NODES];  //substrate stoichiometric values in reverse direction - Sr
int i, j, k;
float maxCncDiff = 0.0f;

printf("\nStoic\n");
//forward and reverse stoichiometric matrices
for(i = 0; i < NUM_NODES * NUM_EDGES; i++)
{
       if(fullStoic[i] < 0)
       {
           fwdStoic[i] = -1 * fullStoic[i];
           revStoic[i] = 0.0f;
       }
       else if(fullStoic[i] > 0)
       {
           fwdStoic[i] = 0.0f;
           revStoic[i] = fullStoic[i];  
       }
       else 
       { 
           fwdStoic[i] = 0.0f;
           revStoic[i] = 0.0f;  
       }
      printf(" \t %f \t %f \n",  fwdStoic[i],revStoic[i]);
}

maxCncDiff = max(&metaCncDiff[0], numMetabolites);

//while(maxCncDiff > tolerance)
//{
for(i = 0; i < numMetabolites ; i++)
   {
       scaleFactor[i] = 1;
}

for(timeStep = 0; timeStep < 1000; timeStep++)
{
// loop through the fluxes and calculate the fwd and reverse flux vectors
    printf("\n fluxes \n");
    for(i = 0; i < numReactions; i++)
    {
        fwdFlux[i] = fwdRate[i];
        revFlux[i] = revRate[i];
        for(j = 0; j < numMetabolites; j++)
        { 
            fwdFlux[i] = fwdFlux[i] * powf(metaCnc[j],fwdStoic[(j * numReactions) + i]); // indices are reversed in matlab code
            revFlux[i] = revFlux[i] * powf(metaCnc[j],revStoic[(j * numReactions) + i]);
        }

    printf("\t %f \t %f \n", fwdFlux[i], revFlux[i]);
    }
    // calculate the metabolic concentration change
    printf("\n meta diff \n");
    for(i = 0; i < numMetabolites ; i++)
    {
        metaCncDiff[i] = 0;
        for(j = 0; j < numReactions; j++)
        {
                metaCncDiff[i] += fullStoic[(i * numReactions) + j] * (fwdFlux[j] - revFlux[j]);
        }
      printf("\t %f ", metaCncDiff[i]);
    }
    printf("\n meta conc \n");
    // Compute new metabolite concentrations
    for(i = 0; i < numMetabolites; i++)
    {
        metaCnc[i] = metaCnc[i] + (scaleFactor[i] * metaCncDiff[i] * t);
    }

    for(i = 0; i < numMetabolites; i++)
    {
        // if concentration is negative, set it to a very small value.
        if(metaCnc[i] < 0)
        {
            metaCnc[i] = .00001;
        }
    printf("\t %f ", metaCnc[i]);
    }

    // Compute new volume scale factor
    printf("\n a \n");
    for(i = 0; i < numMetabolites; i++)
    {
        scaleFactor[i] = epsilon + ( 1 / (abs(metaCncDiff[i]) + epsilon));
        printf("\t %f ", metaCnc[i]);
    }

    //maxCncDiff = max(&metaCncDiff[0], numMetabolites);
}
//} //while


// Calculate total flux
for(i = 0; i < numReactions; i++)
{
   totFlux[i] = fwdFlux[i] - revFlux[i];
}


printf("\n vt \n");
for(i = 0; i < numReactions; i++)
{
   printf("\t %f ", totFlux[i]);
}
printf("\n");
return 0;
}

//Returns the minimum value in the array
float max(float* input, int size)
{
   int i = 0;
   float maxm = abs(*input);

   for(i = 1; i < size; i++)
   {
        maxm = (maxm >= abs(*(input+i))) ? maxm : abs(*(input+i));
   }

return maxm;
}

/*******************************************************************************/

//code to read from the input files
/*if(argc != 0)
{
ReadFile(argv[0], &fwdRate[0]);
ReadFile(argv[1], &revRate[0]);
}*/


/*

// Initialize values
for(i = 0; i < NUM_EDGES; i++)
{

    fwdFlux[i] = 1;
    revFlux[i] = 1;
    //Dummy values - change
    fwdRate[i] = rand() % 10;
    revRate[i] = rand() % 10;
}

for(j = 0; j < NUM_NODES; j++)
{
   scaleFactor[i] = 1;
   metaCnc[j] = 1;
   metaCncDiff[j] = 1;
}
*/



/*
//Outputs
printf("Volume Scale factor \n");

for(i = 0; i < NUM_EDGES; i++)
{
printf("%f \t", scaleFactor[i]);
}
printf("\n");
printf("Forward Flux \n");
for(i = 0; i < NUM_EDGES; i++)
{
printf("%f \t", fwdFlux[i]);
}
printf("\n");
printf("Reverse Flux \n");
for(i = 0; i < NUM_EDGES; i++)
{
printf("%f \t", revFlux[i]);
}
printf("\n");
printf("Metabolic Concetration \n");
for(i = 0; i < NUM_NODES; i++)
{
printf("%f \t", metaCnc[i]);
}
*/


/*****************************************************************/
/*
//Read input files
void ReadFile(char* filename, void* input)
{
    unsigned int i = 0;
    FILE *fp;
    fp = fopen(filename, "r");
    if(fp != NULL)
    {
        while(!feof(fp))
        {
            fscanf(fp, "%f", &input[i]);
            printf("%f \n", input[i]);
            i++;
        }
    }
   fclose(fp);
}
}*/
